library(plotly)
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(datasets)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:plotly':
##
## select
library(glmnet)
## Loading required package: Matrix
## Loaded glmnet 4.1-8
library(ridge)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
Загрузите данные из файла reglab1.txt. Используя функцию lm, постройте регрессию (используйте разные модели). Выберите наиболее подходящую модель, объясните свой выбор.
data <- read.table("reglab1.txt", header = TRUE)
head(data)
## z x y
## 1 2.836772 0.2710104 0.3083308
## 2 4.987167 0.5895982 0.5149134
## 3 6.412325 0.6517442 0.7304530
## 4 4.641998 0.5819827 0.4614002
## 5 2.793941 0.4636884 0.1911021
## 6 2.934520 0.3835542 0.2612945
n <- dim(data)[1]
data_rand <- data[order(runif(n)),]
df_train <- data_rand[1:as.integer(n*0.8),]
df_test <- data_rand[(as.integer(n*0.8)+1):n,]
plot_ly(data = df_train) %>%
add_markers(
z = ~z, x = ~x, y = ~y,
opacity = 0.7,
size = I(100)
) %>%
layout(title = "Training Data")
plot_ly(data = df_test) %>%
add_markers(
z = ~z, x = ~x, y = ~y,
opacity = 0.7,
size = I(100)
) %>%
layout(title = "Test Data")
model1 <- lm(z ~ ., df_train)
predicted_z <- predict(model1, df_test)
model2 <- lm(y ~ ., df_train)
predicted_y <- predict(model2, df_test)
model3 <- lm(x ~ ., df_train)
predicted_x <- predict(model3, df_test)
# Ошибки для моделей
mist_z <- sd(df_test$z - predicted_z)
model1
##
## Call:
## lm(formula = z ~ ., data = df_train)
##
## Coefficients:
## (Intercept) x y
## -0.04412 4.10877 5.01872
cat("Ошибка на тестовых данных: ", mist_z, "\n")
## Ошибка на тестовых данных: 0.3443943
mist_y <- sd(df_test$y - predicted_y)
model2
##
## Call:
## lm(formula = y ~ ., data = df_train)
##
## Coefficients:
## (Intercept) z x
## 0.02944 0.18994 -0.77650
cat("Ошибка на тестовых данных: ", mist_y, "\n")
## Ошибка на тестовых данных: 0.07083249
mist_x <- sd(df_test$x - predicted_x)
model3
##
## Call:
## lm(formula = x ~ ., data = df_train)
##
## Coefficients:
## (Intercept) z y
## 0.04961 0.22330 -1.11505
cat("Ошибка на тестовых данных: ", mist_x, "\n")
## Ошибка на тестовых данных: 0.08171881
Наиболее подходящая модель - модель y, так как меньше ошибка на тестовых данных.
Реализуйте следующий алгоритм для уменьшения количества признаков, используемых для построения регрессии: для каждого k из {0,1,…,d} выбрать подмножество признаков мощности k^1, минимизирующее остаточную сумму квадратов RSS. Используя полученный алгоритм, выберите оптимальное подможество признаков для данных из файла reglab2.txt. Объясните свой выбор. Для генерации всех возможных сочетаний по m элементов из некоторого множества x можно использовать функцию combn(x, m, …).
data <- read.table("reglab2.txt", header = TRUE)
head(data)
## y x1 x2 x3 x4
## 1 3.4697200 0.23362753 0.83554899 0.1029653 0.45742768
## 2 0.7684485 0.11791994 0.09054379 0.2587784 0.28395097
## 3 2.8803737 0.09152018 0.79759245 0.1985277 0.69928681
## 4 3.7454851 0.87672239 0.06293480 0.6154150 0.17605768
## 5 1.8539662 0.20740592 0.30348994 0.7759667 0.66735068
## 6 1.2952686 0.04407584 0.36008230 0.6572325 0.04344568
compute_RSS <- function(features, data) {
formula <- as.formula(paste("y ~", paste(features, collapse = " + ")))
model <- lm(formula, data = data)
RSS <- sum(model$residuals^2)
return(RSS)
}
find_optimal_subset_lm <- function(data) {
p <- ncol(data) - 1
features <- colnames(data)[-1]
best_RSS <- Inf
best_subset <- NULL
all_results <- list()
# Итерация по всем возможным размерам подмножества
for (m in 1:p) {
subsets <- combn(features, m, simplify = TRUE)
# Итерация по всем возможным сочетаниям признаков данного размера
for (i in 1:ncol(subsets)) {
current_subset <- subsets[, i]
current_RSS <- compute_RSS(current_subset, data)
all_results[[length(all_results) + 1]] <- list(subset = current_subset, RSS = current_RSS)
if (current_RSS < best_RSS) {
best_RSS <- current_RSS
best_subset <- current_subset
}
}
}
return(list(all_results = all_results, optimal_subset = best_subset, optimal_RSS = best_RSS))
}
result <- find_optimal_subset_lm(data)
cat("Все подмножества и их RSS:\n")
## Все подмножества и их RSS:
for (i in seq_along(result$all_results)) {
cat("Подмножество:", paste(result$all_results[[i]]$subset, collapse = ", "), " | RSS:", result$all_results[[i]]$RSS, "\n")
}
## Подмножество: x1 | RSS: 157.2198
## Подмножество: x2 | RSS: 268.2458
## Подмножество: x3 | RSS: 393.4905
## Подмножество: x4 | RSS: 394.5905
## Подмножество: x1, x2 | RSS: 0.5379617
## Подмножество: x1, x3 | RSS: 156.3541
## Подмножество: x1, x4 | RSS: 157.2193
## Подмножество: x2, x3 | RSS: 267.7955
## Подмножество: x2, x4 | RSS: 267.8061
## Подмножество: x3, x4 | RSS: 393.4587
## Подмножество: x1, x2, x3 | RSS: 0.3322662
## Подмножество: x1, x2, x4 | RSS: 0.3619682
## Подмножество: x1, x3, x4 | RSS: 156.3483
## Подмножество: x2, x3, x4 | RSS: 267.4415
## Подмножество: x1, x2, x3, x4 | RSS: 0.1928635
cat("\nОптимальное подмножество признаков:", paste(result$optimal_subset, collapse = ", "), "\n")
##
## Оптимальное подмножество признаков: x1, x2, x3, x4
cat("Максимальная остаточная сумма квадратов (RSS):", result$optimal_RSS, "\n")
## Максимальная остаточная сумма квадратов (RSS): 0.1928635
Загрузите данные из файла cygage.txt. Постройте регрессию, выражающую зависимость возраста исследуемых отложений от глубины залегания, используя веса наблюдений. Оцените качество построенной модели.
data <- read.table("cygage.txt", header = TRUE)
head(data)
## calAge Depth Weight
## 1 0 0 1.0
## 2 3707 77 0.1
## 3 4150 90 0.1
## 4 5350 107 0.1
## 5 4500 168 0.1
## 6 7260 217 0.1
plot(data, col = 'blue')
model <- lm(calAge ~ Depth, data, weights = data$Weight)
summary <- summary(model)
mid_error = mean(summary$residuals^2)
cat("Среднеквадратичная ошибка (MSE):", mid_error, "\n")
## Среднеквадратичная ошибка (MSE): 227161.2
Загрузите данные Longley (макроэкономические данные). Данные состоят из 7 экономических переменных, наблюдаемых с 1947 по 1962 годы (n=16): GNP.deflator - дефлятор цен, GNP - валовой национальный продукт, Unemployed – число безработных Armed.Forces – число людей в армии Population – население, возраст которого старше 14 лет Year - год Employed – количество занятых Построить регрессию lm(Employed ~ .). Исключите из набора данных longley переменную “Population”. Разделите данные на тестовую и обучающую выборки равных размеров случайным образом. Постройте гребневую регрессию для значений lambda=10^(-3+0.2*i), i=0,…,25, подсчитайте ошибку на тестовой и обучающей выборке для данных значений λ, постройте графики. Объясните полученные результаты.
data(longley)
head(longley)
## GNP.deflator GNP Unemployed Armed.Forces Population Year Employed
## 1947 83.0 234.289 235.6 159.0 107.608 1947 60.323
## 1948 88.5 259.426 232.5 145.6 108.632 1948 61.122
## 1949 88.2 258.054 368.2 161.6 109.773 1949 60.171
## 1950 89.5 284.599 335.1 165.0 110.929 1950 61.187
## 1951 96.2 328.975 209.9 309.9 112.075 1951 63.221
## 1952 98.1 346.999 193.2 359.4 113.270 1952 63.639
model <- lm(Employed ~ ., longley)
summary <- summary(model)
mid_error = mean(summary$residuals^2)
cat("Среднеквадратичная ошибка (MSE):", mid_error, "\n")
## Среднеквадратичная ошибка (MSE): 0.0522765
longley$Population <- NULL
n <- dim(longley)[1]
longley_rand <- longley[order(runif(n)),]
df_train <- longley_rand[1:as.integer(n*0.8),]
df_test <- longley_rand[(as.integer(n*0.8)+1):n,]
lambda_values <- 10^seq(-3, 2, by = 0.2)
train_error <- numeric(length(lambda_values))
test_error <- numeric(length(lambda_values))
for (i in seq_along(lambda_values)) {
lambda <- lambda_values[i]
model <- lm.ridge(Employed ~ ., df_train, lambda = lambda)
coef_values <- as.matrix(coef(model))
train_pred <- as.matrix(cbind(1, df_train[, -ncol(df_train)])) %*% coef_values
test_pred <- as.matrix(cbind(1, df_test[, -ncol(df_test)])) %*% coef_values
train_error[i] <- mean((train_pred - df_train$Employed)^2)
test_error[i] <- mean((test_pred - df_test$Employed)^2)
}
plot(log10(lambda_values), train_error, type = "l", col = "blue", xlab = "log10(lambda)", ylab = "MSE", main = "Гребневая регрессия")
lines(log10(lambda_values), test_error, type = "l", col = "red")
legend("topright", legend = c("Обучающая выборка", "Тестовая выборка"), col = c("blue", "red"), lty = 1)
Из графика видно, что среднеквадратичная ошибка на обучающей выборке больше, чем на тестовой. Причем при увеличении lambda, ошибка растет.
Загрузите данные EuStockMarkets из пакета « datasets». Данные содержат ежедневные котировки на момент закрытия фондовых бирж: Germany DAX (Ibis), Switzerland SMI, France CAC, и UK FTSE. Постройте на одном графике все кривые изменения котировок во времени. Постройте линейную регрессию для каждой модели в отдельности и для всех моделей вместе. Оцените, какая из бирж имеет наибольшую динамику.
data(EuStockMarkets)
head(EuStockMarkets)
## DAX SMI CAC FTSE
## [1,] 1628.75 1678.1 1772.8 2443.6
## [2,] 1613.63 1688.5 1750.5 2460.2
## [3,] 1606.51 1678.6 1718.0 2448.2
## [4,] 1621.04 1684.1 1708.1 2470.4
## [5,] 1618.16 1686.6 1723.1 2484.7
## [6,] 1610.61 1671.6 1714.3 2466.8
plot(EuStockMarkets[,1], type = "l", xlab="Time", ylab="Closing Price", col="black")
lines(EuStockMarkets[,2], type = "l", col="red")
lines(EuStockMarkets[,3], type = "l", col="green")
lines(EuStockMarkets[,4], type = "l", col="blue")
legend("topleft", legend = colnames(EuStockMarkets), col = 1:4, lty = 1)
model1 <- lm(DAX ~ time(EuStockMarkets), EuStockMarkets)
model2 <- lm(SMI ~ time(EuStockMarkets), EuStockMarkets)
model3 <- lm(CAC ~ time(EuStockMarkets), EuStockMarkets)
model4 <- lm(FTSE ~ time(EuStockMarkets), EuStockMarkets)
model5 <- lm(DAX+SMI+CAC+FTSE ~ time(EuStockMarkets), EuStockMarkets)
cat("-----------Регрессия для DAX----------- \n")
## -----------Регрессия для DAX-----------
model1
##
## Call:
## lm(formula = DAX ~ time(EuStockMarkets), data = EuStockMarkets)
##
## Coefficients:
## (Intercept) time(EuStockMarkets)
## -894557.9 449.7
cat("-----------Регрессия для SMI----------- \n")
## -----------Регрессия для SMI-----------
model2
##
## Call:
## lm(formula = SMI ~ time(EuStockMarkets), data = EuStockMarkets)
##
## Coefficients:
## (Intercept) time(EuStockMarkets)
## -1428160.2 717.5
cat("-----------Регрессия для CAC----------- \n")
## -----------Регрессия для CAC-----------
model3
##
## Call:
## lm(formula = CAC ~ time(EuStockMarkets), data = EuStockMarkets)
##
## Coefficients:
## (Intercept) time(EuStockMarkets)
## -405915.3 204.6
cat("-----------Регрессия для FTSE----------- \n")
## -----------Регрессия для FTSE-----------
model4
##
## Call:
## lm(formula = FTSE ~ time(EuStockMarkets), data = EuStockMarkets)
##
## Coefficients:
## (Intercept) time(EuStockMarkets)
## -865200.4 435.5
cat("-----------Регрессия для всех----------- \n")
## -----------Регрессия для всех-----------
model5
##
## Call:
## lm(formula = DAX + SMI + CAC + FTSE ~ time(EuStockMarkets), data = EuStockMarkets)
##
## Coefficients:
## (Intercept) time(EuStockMarkets)
## -3593834 1807
Из результатов видим, что наибольшую динамику имеет SMI, так как у неё самый большой коэффициент при параметре time.
Загрузите данные JohnsonJohnson из пакета «datasets». Данные содержат поквартальную прибыль компании Johnson & Johnson с 1960 по 1980 гг. Постройте на одном графике все кривые изменения прибыли во времени. Постройте линейную регрессию для каждого квартала в отдельности и для всех кварталов вместе. Оцените, в каком квартале компания имеет наибольшую и наименьшую динамику доходности. Сделайте прогноз по прибыли в 2016 году во всех кварталах и в среднем по году.
data(JohnsonJohnson)
JohnsonJohnson
## Qtr1 Qtr2 Qtr3 Qtr4
## 1960 0.71 0.63 0.85 0.44
## 1961 0.61 0.69 0.92 0.55
## 1962 0.72 0.77 0.92 0.60
## 1963 0.83 0.80 1.00 0.77
## 1964 0.92 1.00 1.24 1.00
## 1965 1.16 1.30 1.45 1.25
## 1966 1.26 1.38 1.86 1.56
## 1967 1.53 1.59 1.83 1.86
## 1968 1.53 2.07 2.34 2.25
## 1969 2.16 2.43 2.70 2.25
## 1970 2.79 3.42 3.69 3.60
## 1971 3.60 4.32 4.32 4.05
## 1972 4.86 5.04 5.04 4.41
## 1973 5.58 5.85 6.57 5.31
## 1974 6.03 6.39 6.93 5.85
## 1975 6.93 7.74 7.83 6.12
## 1976 7.74 8.91 8.28 6.84
## 1977 9.54 10.26 9.54 8.73
## 1978 11.88 12.06 12.15 8.91
## 1979 14.04 12.96 14.85 9.99
## 1980 16.20 14.67 16.02 11.61
time_labels <- time(JohnsonJohnson)[seq(from = 1, to = length(JohnsonJohnson), by = 4)]
col_names <- c('Qtr1','Qtr2','Qtr3','Qtr4')
Qtr1 <- JohnsonJohnson[seq(from = 1, to = length(JohnsonJohnson), by = 4)]
Qtr2 <- JohnsonJohnson[seq(from = 2, to = length(JohnsonJohnson), by = 4)]
Qtr3 <- JohnsonJohnson[seq(from = 3, to = length(JohnsonJohnson), by = 4)]
Qtr4 <- JohnsonJohnson[seq(from = 4, to = length(JohnsonJohnson), by = 4)]
plot(time_labels, Qtr1, type = "l", xlab = "Время", ylab = "Кривая", col = "black")
lines(time_labels, Qtr2, type = "l", col = "red")
lines(time_labels, Qtr3, type = "l", col = "green")
lines(time_labels, Qtr4, type = "l", col = "blue")
legend("topleft", legend = col_names, col = 1:4, lty = 1)
model1 <- lm(Qtr1 ~ time_labels, JohnsonJohnson)
model2 <- lm(Qtr2 ~ time_labels, JohnsonJohnson)
model3 <- lm(Qtr3 ~ time_labels, JohnsonJohnson)
model4 <- lm(Qtr4 ~ time_labels, JohnsonJohnson)
model5 <- lm(Qtr1 + Qtr2 + Qtr3 + Qtr4 ~ time_labels, JohnsonJohnson)
cat("-----------Регрессия для Qtr1----------- \n")
## -----------Регрессия для Qtr1-----------
model1
##
## Call:
## lm(formula = Qtr1 ~ time_labels, data = JohnsonJohnson)
##
## Coefficients:
## (Intercept) time_labels
## -1364.282 0.695
cat("-----------Регрессия для Qtr2----------- \n")
## -----------Регрессия для Qtr2-----------
model2
##
## Call:
## lm(formula = Qtr2 ~ time_labels, data = JohnsonJohnson)
##
## Coefficients:
## (Intercept) time_labels
## -1345.0727 0.6853
cat("-----------Регрессия для Qtr3----------- \n")
## -----------Регрессия для Qtr3-----------
model3
##
## Call:
## lm(formula = Qtr3 ~ time_labels, data = JohnsonJohnson)
##
## Coefficients:
## (Intercept) time_labels
## -1382.3170 0.7044
cat("-----------Регрессия для Qtr4----------- \n")
## -----------Регрессия для Qtr4-----------
model4
##
## Call:
## lm(formula = Qtr4 ~ time_labels, data = JohnsonJohnson)
##
## Coefficients:
## (Intercept) time_labels
## -1049.5828 0.5349
cat("-----------Регрессия для всех----------- \n")
## -----------Регрессия для всех-----------
model5
##
## Call:
## lm(formula = Qtr1 + Qtr2 + Qtr3 + Qtr4 ~ time_labels, data = JohnsonJohnson)
##
## Coefficients:
## (Intercept) time_labels
## -5141.25 2.62
Из результатов видим, что наибольшая динамика наблюдается в 3 квартале, так как у него самый большой коэффициент при параметре time. Наименьшая же динамика наблюдается в 4 квартале.
predict1 = coef(model1)[1]+coef(model1)[2]*2016
predict2 = coef(model2)[1]+coef(model2)[2]*2016
predict3 = coef(model3)[1]+coef(model3)[2]*2016
predict4 = coef(model4)[1]+coef(model4)[2]*2016
predict5 = (coef(model5)[1]+coef(model5)[2]*2016)/4
cat("Прогноз прибыли в 2016 году для Qtr1: ", predict1, "\n")
## Прогноз прибыли в 2016 году для Qtr1: 36.75964
cat("Прогноз прибыли в 2016 году для Qtr2: ", predict2, "\n")
## Прогноз прибыли в 2016 году для Qtr2: 36.48945
cat("Прогноз прибыли в 2016 году для Qtr3: ", predict3, "\n")
## Прогноз прибыли в 2016 году для Qtr3: 37.65394
cat("Прогноз прибыли в 2016 году для Qtr4: ", predict4, "\n")
## Прогноз прибыли в 2016 году для Qtr4: 28.79391
cat("Прогноз прибыли в 2016 году в среднем по году: ", predict5, "\n")
## Прогноз прибыли в 2016 году в среднем по году: 34.92424
Загрузите данные sunspot.year из пакета «datasets». Данные содержат количество солнечных пятен с 1700 по 1988 гг. Постройте на графике кривую изменения числа солнечных пятен во времени. Постройте линейную регрессию для данных.
data(sunspot.year)
sunspot.year
## Time Series:
## Start = 1700
## End = 1988
## Frequency = 1
## [1] 5.0 11.0 16.0 23.0 36.0 58.0 29.0 20.0 10.0 8.0 3.0 0.0
## [13] 0.0 2.0 11.0 27.0 47.0 63.0 60.0 39.0 28.0 26.0 22.0 11.0
## [25] 21.0 40.0 78.0 122.0 103.0 73.0 47.0 35.0 11.0 5.0 16.0 34.0
## [37] 70.0 81.0 111.0 101.0 73.0 40.0 20.0 16.0 5.0 11.0 22.0 40.0
## [49] 60.0 80.9 83.4 47.7 47.8 30.7 12.2 9.6 10.2 32.4 47.6 54.0
## [61] 62.9 85.9 61.2 45.1 36.4 20.9 11.4 37.8 69.8 106.1 100.8 81.6
## [73] 66.5 34.8 30.6 7.0 19.8 92.5 154.4 125.9 84.8 68.1 38.5 22.8
## [85] 10.2 24.1 82.9 132.0 130.9 118.1 89.9 66.6 60.0 46.9 41.0 21.3
## [97] 16.0 6.4 4.1 6.8 14.5 34.0 45.0 43.1 47.5 42.2 28.1 10.1
## [109] 8.1 2.5 0.0 1.4 5.0 12.2 13.9 35.4 45.8 41.1 30.1 23.9
## [121] 15.6 6.6 4.0 1.8 8.5 16.6 36.3 49.6 64.2 67.0 70.9 47.8
## [133] 27.5 8.5 13.2 56.9 121.5 138.3 103.2 85.7 64.6 36.7 24.2 10.7
## [145] 15.0 40.1 61.5 98.5 124.7 96.3 66.6 64.5 54.1 39.0 20.6 6.7
## [157] 4.3 22.7 54.8 93.8 95.8 77.2 59.1 44.0 47.0 30.5 16.3 7.3
## [169] 37.6 74.0 139.0 111.2 101.6 66.2 44.7 17.0 11.3 12.4 3.4 6.0
## [181] 32.3 54.3 59.7 63.7 63.5 52.2 25.4 13.1 6.8 6.3 7.1 35.6
## [193] 73.0 85.1 78.0 64.0 41.8 26.2 26.7 12.1 9.5 2.7 5.0 24.4
## [205] 42.0 63.5 53.8 62.0 48.5 43.9 18.6 5.7 3.6 1.4 9.6 47.4
## [217] 57.1 103.9 80.6 63.6 37.6 26.1 14.2 5.8 16.7 44.3 63.9 69.0
## [229] 77.8 64.9 35.7 21.2 11.1 5.7 8.7 36.1 79.7 114.4 109.6 88.8
## [241] 67.8 47.5 30.6 16.3 9.6 33.2 92.6 151.6 136.3 134.7 83.9 69.4
## [253] 31.5 13.9 4.4 38.0 141.7 190.2 184.8 159.0 112.3 53.9 37.5 27.9
## [265] 10.2 15.1 47.0 93.8 105.9 105.5 104.5 66.6 68.9 38.0 34.5 15.5
## [277] 12.6 27.5 92.5 155.4 154.7 140.5 115.9 66.6 45.9 17.9 13.4 29.2
## [289] 100.2
plot(sunspot.year, col='blue')
model <- lm(sunspot.year[seq(1,length(sunspot.year),1)] ~ seq(1700,1988,1), sunspot.year)
cat("-----------Регрессия----------- \n")
## -----------Регрессия-----------
model
##
## Call:
## lm(formula = sunspot.year[seq(1, length(sunspot.year), 1)] ~
## seq(1700, 1988, 1), data = sunspot.year)
##
## Coefficients:
## (Intercept) seq(1700, 1988, 1)
## -130.42119 0.09709
summary <- summary(model)
mid_error = mean(summary$residuals^2)
cat("Среднеквадратичная ошибка (MSE):", mid_error, "\n")
## Среднеквадратичная ошибка (MSE): 1487.204
Загрузите данные из файла пакета «UKgas.scv». Данные содержат объемы ежеквартально потребляемого газа в Великобритании с 1960 по 1986 гг. Постройте линейную регрессию для каждого квартала в отдельности и для всех кварталов вместе. Оцените, в каком квартале потребление газа имеет наибольшую и наименьшую динамику доходности. Сделайте прогноз по потреблению газа в 2016 году во всех кварталах и в среднем по году.
data <- read.csv("UKgas.csv", stringsAsFactors = TRUE)
head(data)
## X time UKgas
## 1 1 1960.00 160.1
## 2 2 1960.25 129.7
## 3 3 1960.50 84.8
## 4 4 1960.75 120.1
## 5 5 1961.00 160.1
## 6 6 1961.25 124.9
Qtr1 <- UKgas[seq(from = 1, to = dim(data)[1], by = 4)]
Qtr2 <- UKgas[seq(from = 2, to = dim(data)[1], by = 4)]
Qtr3 <- UKgas[seq(from = 3, to = dim(data)[1], by = 4)]
Qtr4 <- UKgas[seq(from = 4, to = dim(data)[1], by = 4)]
model1 <- lm(Qtr1 ~ time[seq(from = 1, to = dim(data)[1], by = 4)], data)
model2 <- lm(Qtr2 ~ time[seq(from = 2, to = dim(data)[1], by = 4)], data)
model3 <- lm(Qtr3 ~ time[seq(from = 3, to = dim(data)[1], by = 4)], data)
model4 <- lm(Qtr4 ~ time[seq(from = 4, to = dim(data)[1], by = 4)], data)
model5 <- lm(Qtr1 + Qtr2 + Qtr3 + Qtr4 ~ time[seq(from = 1, to = dim(data)[1], by = 4)], data)
cat("-----------Регрессия для Qtr1----------- \n")
## -----------Регрессия для Qtr1-----------
model1
##
## Call:
## lm(formula = Qtr1 ~ time[seq(from = 1, to = dim(data)[1], by = 4)],
## data = data)
##
## Coefficients:
## (Intercept)
## -78854.23
## time[seq(from = 1, to = dim(data)[1], by = 4)]
## 40.22
cat("-----------Регрессия для Qtr2----------- \n")
## -----------Регрессия для Qtr2-----------
model2
##
## Call:
## lm(formula = Qtr2 ~ time[seq(from = 2, to = dim(data)[1], by = 4)],
## data = data)
##
## Coefficients:
## (Intercept)
## -35297.23
## time[seq(from = 2, to = dim(data)[1], by = 4)]
## 18.04
cat("-----------Регрессия для Qtr3----------- \n")
## -----------Регрессия для Qtr3-----------
model3
##
## Call:
## lm(formula = Qtr3 ~ time[seq(from = 3, to = dim(data)[1], by = 4)],
## data = data)
##
## Coefficients:
## (Intercept)
## -15403.73
## time[seq(from = 3, to = dim(data)[1], by = 4)]
## 7.89
cat("-----------Регрессия для Qtr4----------- \n")
## -----------Регрессия для Qtr4-----------
model4
##
## Call:
## lm(formula = Qtr4 ~ time[seq(from = 4, to = dim(data)[1], by = 4)],
## data = data)
##
## Coefficients:
## (Intercept)
## -59112.72
## time[seq(from = 4, to = dim(data)[1], by = 4)]
## 30.14
cat("-----------Регрессия для всех----------- \n")
## -----------Регрессия для всех-----------
model5
##
## Call:
## lm(formula = Qtr1 + Qtr2 + Qtr3 + Qtr4 ~ time[seq(from = 1, to = dim(data)[1],
## by = 4)], data = data)
##
## Coefficients:
## (Intercept)
## -188636.85
## time[seq(from = 1, to = dim(data)[1], by = 4)]
## 96.29
Из результатов видим, что наибольшая динамика наблюдается в 1 квартале, так как у него самый большой коэффициент при параметре time. Наименьшая же динамика наблюдается в 3 квартале.
predict1 = coef(model1)[1]+coef(model1)[2]*2016
predict2 = coef(model2)[1]+coef(model2)[2]*2016
predict3 = coef(model3)[1]+coef(model3)[2]*2016
predict4 = coef(model4)[1]+coef(model4)[2]*2016
predict5 = (coef(model5)[1]+coef(model5)[2]*2016)/4
cat("Прогноз прибыли в 2016 году для Qtr1: ", predict1, "\n")
## Прогноз прибыли в 2016 году для Qtr1: 2230.936
cat("Прогноз прибыли в 2016 году для Qtr2: ", predict2, "\n")
## Прогноз прибыли в 2016 году для Qtr2: 1072.375
cat("Прогноз прибыли в 2016 году для Qtr3: ", predict3, "\n")
## Прогноз прибыли в 2016 году для Qtr3: 501.9919
cat("Прогноз прибыли в 2016 году для Qtr4: ", predict4, "\n")
## Прогноз прибыли в 2016 году для Qtr4: 1654.785
cat("Прогноз прибыли в 2016 году в среднем по году: ", predict5, "\n")
## Прогноз прибыли в 2016 году в среднем по году: 1372.787
Загрузите данные cars из пакета «datasets». Данные содержат зависимости тормозного пути автомобиля (футы) от его скорости (мили в час). Данные получены в 1920 г. Постройте регрессионную модель и оцените длину тормозного пути при скорости 40 миль в час.
data(cars)
head(cars)
## speed dist
## 1 4 2
## 2 4 10
## 3 7 4
## 4 7 22
## 5 8 16
## 6 9 10
model <- lm(dist ~ speed, cars)
cat("-----------Регрессия----------- \n")
## -----------Регрессия-----------
model
##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Coefficients:
## (Intercept) speed
## -17.579 3.932
predict = coef(model)[1]+coef(model)[2]*40
cat("Прогноз длины тормозного пути при скорости 40 миль в час: ", predict, "\n")
## Прогноз длины тормозного пути при скорости 40 миль в час: 139.7173